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1,  Introduction 


In  the  area  of  wave  propagation  in  random  media,  most  workers  assumed 
that  the  fluctuations  are  independent  of  time.  For  a  short-range  propa¬ 
gation  at  high  frequency,  the  omission  of  temporal  fluctuation  may  not 
cause  serious  errors.  However,  in  the  long-range  transmission  or  in  the 
propagation  of  pulses,  the  temporal  fluctuations  of  media  become  important 
and  should  be  taken  into  account.  For  physical  reason  a  time-dependent 
random  medium  will  be  called  a  turbulent  medium.  For  earlier  papers  of 
basic  importance,  one  is  referred  to  the  classics  by  Keller  [1,  2]. 

In  this  paper  we  shall  develop  a  general  effective  method  of  computing 
the  statistics  of  a  pulsed  wave  travelling  in  a  turbulent  medium.  To  any 
experienced  worker  in  this  field,  it  seems  clear  that  the  goal  is  unattainable 
in  the  strong  fluctuation  case.  Therefore  our  study  will  be  confined  to 
the  case  of  weak  fluctuation.  Then  it  is  possible  to  devise  a  method  which 
will  account  for  the  accumulative  effect  of  weak  fluctuations  in  the  long 
run.  If  the  governing  equations  were  ordinary,  instead  of  partial  differen¬ 
tial  equations,  the  powerful  method  of  diffusion  approximation  developed 
by  Khasiminskii [3] ,  and  Papanicolaou  [4]  would  have  been  an  answer.  Here, 
though  not  directly  applicable,  it  will  also  play  an  important  role  as  to 
be  seen.  Our  method  consists  of  three  basic  steps:  the  deterministic 
method  of  characteristics  for  PDE’s  (partial  differential  equations),  the 
diffusion  approximation  for  stochastic  ODE's  (ordinary  differential  equations), 
and  the  method  of  functional  integrals  [5] ,  in  that  order.  As  the  first 
step,  we  seek  a  progressing  wave  solution  which  enables  us  to  reduce  a 
stochastic  PDE  to  a  system  of  stochastic  ODE's,  to  which  the  diffusion 
approximation  is  applicable.  By  invoking  the  diffusion  limit  theorems 
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mentioned  above,  the  amplitude  and  phase  fluctuations  may  be  described 
asymptotically  in  terms  of  some  known  diffusion  processes.  Thereby  the 
statistics  of  the  original  random  wave  function  can  be  computed  approxi¬ 
mately  by  evaluating  some  average  functionals  of  the  known  diffusion  pro¬ 
cesses.  This  last  step  may  be  viewed  as  the  evaluation  of  certain 
functional  integrals  over  the  probability  distributions  of  the  diffusion 
processes.  In  probability  theory,  the  last  two  steps  correspond  to  the 
celebrated  invariance  principle  of  Ponsker  [6],  in  a  generalized  situation. 
We  wish  to  point  out  that,  even  after  an  asymptotic  reduction,  the  evalua¬ 
tion  of  a  functional  integral  is  by  no  means  easy.  But  the  reduced  problem 
is  much  simpler  than  the  original  one  which  is  totally  intractable.  The 
proposed  method  is  well  suited  for  treating  the  pulse  propagation  problem 
governed  by  a  random  hyperbolic  system.  It  yields  a  statistical  description 
of  the  physically  meaningful  quantities,  such  as  the  ray,  amplitude  and 
phase,  in  terms  of  which  the  statistical  law  for  the  wave  function  may  be 
determined.  This  allows  us,  in  principle,  to  compute  more  complicated 
statistics,  other  than  the  moments  of  the  solution  as  commonly  done. 

The  text  of  the  paper  consists  of  three  sections.  In  Section  2  we 
develop  the  basic  ideas  via  a  simple  model  equation.  The  first-order 
random  PDF,  is  to  be  solved  by  the  method  of  characteristics.  By  applying 
the  appropriate  limit  theorems,  the  fluctuations  of  amplitude  and  phase 
functions  may  be  determined  along  the  characteristic  curves  in  an  asymptotic 
sense.  As  examples,  some  statistical  functionals  are  evaluated.  The 
general  method  of  progressing  waves  for  a  random  hyperbolic  system  is 
presented  in  Section  3.  Due  to  a  marked  difference  between  one  and  higher 
dimensional  cases,  they  will  be  discussed  separately.  There  the  three 
basic  steps  stated  above  will  be  carried  out  in  detail.  To  illustrate  the 


3 


computational  procedures,  the  method  is  applied  to  problems  in  the  acoustic 
wave  propagation  through  a  turbulent  medium.  This  is  done  in  Section  4. 

To  avoid  obscuring  the  basic  ideas  involved,  some  technical  details  con¬ 
cerning  two  main  theorems  that  we  rely  on  heavily  are  omitted  in  the  text, 
and  they  are  presented  in  the  appendix. 
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2.  Simple  Model  Equation 

To  illustrate  the  basic  ideas,  let  us  consider  a  first-order  partial 
differential  equation  with  a  random  coefficient. 

[3  +  n£(t,x,u>)3  ]u  =  b(t ,x)u,  t  >  0, 

(2.1)  r  t  X 

\  u(0,x,u))  =  a(x)  ,  x  e  R  . 

Here  r)C>  for  each  e>0,  is  a  known  random  function  (field)  of  teE*’  and  xeF  , 
defined  over  a  probability  space  (ft,F,P)  with  the  sample  space  ft  containing 
the  sample  point  u>.  The  given  functions  a  :  1R  ]R  and  b:  R+  *  R  +  IR  or  I 
are  deterministic  and  smooth.  The  solution  u(t,x,u))  is  thus  a  random  field. 
The  statistical  law  and  averages  of  u  are  of  our  primary  concern.  By  the 

standard  notation,  a  statistical  average  of  {•}  will  be  denoted  by  E{-}. 

3  3 

Also,  in  this  paper,  3  ,  3  denote  the  partial  derivatives  ttt  ,  ,  etc. 

t  X  dt  dX 

As  is  well-known,  the  system  (2.1)  can  be  integrated  by  introducing 
the  characteristic  equation: 

(  T7  =  ne(s>y>w),  o<s<t, 

(2.2) 

^  y (t)  =  x, 

which  is  a  stochastic  system. 

Let 

(2.3)  y  =  £^(t,x,co)  ,  o<s<t 

be  a  random  solution  of  (2.2),  which  defines  a  characteristic  curve  ^(co) 
passing  through  (t,x).  Along  such  a  curve,  the  system  (2.1)  can  be  integrated 
to  give  the  "wave  function": 

£  £ 

(2.4)  u(t,x,oo)  =  a[£g(t,x,u))  ]  exp{0  (t,x,oo)}, 


where 
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(2.5)  0e(t,x,u))  =  /Jjb[s,5g(t,x,u))]ds, 

will  be  called  a  "phase  function",  (which  may  not  be  physically  meaningful 
as  shown  in  the  next  section).  We  note  that  the  representation  of  solution 

0 

(2.4)  is  not  constructive,  since,  given  n  ,  the  probability  distribution  of 

0 

can  not  be  determined  easily  in  case  of  strong  fluctuations.  Suppose  that 

(2.6)  ne(t,x,w)  =  n  +  en(t,x,o)), 

where  fj  is  a  constant  and  0<e<<l  is  the  small  scale  of  fluctuations.  Now 
we  let 


(2.7) 


y(s)  =  £s(t,x)  +  z(s) , 
£s(t,x)  =  x  -  n(t-s) . 


Then,  in  view  of  (2.2),  zg  =  z(t-s)  satisfies 


dz 


(2.8) 


dT=  Zt(s*zs*u)’ 


zo  x* 


where 


(2.9)  Zt(s,zs,u))  =  n(t-s,£;t_s+zs,co) . 

To  the  system  (2.8),  we  can  apply  Theorem  A.l  to  get  a  diffusion  approxi¬ 
mation.  To  this  end  we  assume  that  the  random  field  n  is  homogeneous  and 
stationary  such  that 


En  =  0  , 

(2.10) 

En(t,x)n(s,y)  =  Y (t-s ,x-y) , 
and  the  following  limits  exist 

3  =  V2  =  lim  i  fj.  /T  Y[t-s,n(t-s)]  dsdt , 

T  -KO  1  U  L 

(2.H)  i  1  T  fT 

K  =  lim  1  flQ  /*  Yx[t-s,n(t-s)]  dsdt  , 

T-*» 

with  Yx(t,x)  =  ||  Y(t,x) . 
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Then,  under  a  strong  mixing  condition  among  others,  one  can  apply  the 

Khasiminskii  type  of  limit  theorem  to  the  system  (2.8),  as  shown  in  the 

appendix.  According  to  Theorem  A.l  there,  we  can  assert  that,  as  e  4  0, 

2  2 

t>s  t  o°  with  x  =  e  t  and  a  =  e  t  held  fixed,  we  have 
(2 • 12)  =  Cg(t,x)  =>  =  ^(x),  0<o<x 

where  {^o,a<x}  is  a  Wiener  process  with  the  mean  and  variance  parameters  k 
and  v >  respectively,  and  =  0.  Let  (w(t),  t>0}  be  the  standard  Wiener 

process  or  the  Brownian  motion  with  w(0)  =  o.  Then  ^  will  satisfy  the 
ltd  equation 

d£  =  Kda  +  vdw(x-a) ,  a<x  , 

U.13)  ° 

^  =  0,  with  v  =8* 

Therefore,  in  the  diffusion  limit,  the  solution  of  (2.1),  noting  (2.3),  can 
be  approximated  asymptotically  as  follows: 

(2.14)  y(s)  =  £g(t,x)  ~  £s(t,x), 
where,  in  view  of  (2.7)  and  (2.12), 

(2.15)  £g(t,x)  =  x  -  (c+e^c)(t-s)  +  evw(t-s). 

This  shows  clearly  that  the  original  characteristic 

(2.16)  if  :  y  =  d(t,x),  0<s<t , 

l  j  A  O 

is  close  to  the  asymptotic  characteristic 

(2.17)  Tt>x  :  y  =  Cg(t,x)  ,  0<s<t  , 

as  given  explicitly  by  (2.15).  It  contains  a  known  fluctuation,  a  Brownian 
motion  with  the  parameter  v  about  an  effective  (in  general,  not  the  un- 
pertubed)  characteristic 

(2.18)  rt  x  :  y  =  ?s(t,x)  =  x  -  (fj+e2|< )  (t-s) . 
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Suppose  that  the  functions  a  and  b  in  the  system  (2.1)  are  so  smooth 
as  required  by  Theorem  A. 2.  Then  this  theorem  implies  that  the  random 
solution  (2.4)  yields 

u(t,x)  ~  a[£  (t,x)J  exp{0 (t ,x) }  , 

(2.19) 

9(t,x)  =  /Q  b[s,f;s(t,x)]ds. 

Here,  for  brevity,  the  sample  point  u)  in  the  argument  of  all  the  relevent 
functions  is  surpressed.  This  will  often  be  done  hereafter  so  long  as  there 
is  no  confusion.  In  view  of  (2.15),  the  result  (2.19)  shows  that  the  wave 
function  is  in  fact  a  functional  of  Wiener  process  or  Wiener  functional  in 
short.  Hence  the  computation  of  asymptotic  statistics  of  the  wave  function 
u  amounts  to  the  evaluation  of  certain  functional  integrals  with  respect  to 
the  Wiener  measure.  Furthermore  the  statistical  laws  of  the  amplitude  and 
phase  become  completely  known. 

As  an  example,  let  b  be  imaginary, 

(2.20)  b(t,x)  =  i  <f>(t,x),  4>  real. 

Let  <u>  =  Eu  denote  the  mean  or  coherent  wave  function.  Then 

<u>  ~  <v>  =  Ew(a[^(t,x)  -  evw(t)] 

^2'21^  x  exp  i  /q  <j>[s,£(t-s,x)  -  evw(t-s)]ds)  , 

where  E^  denotes  the  expectation  with  respect  to  the  Wiener  process  w(t). 

We  note  that  the  average  Wiener  functional  in  (2.21)  is  just  a  modified 

Feynman-Kac  formula  along  the  effective  characteristic  T  .  It  can  be 

t ,  x 

shown,  similar  to  the  conventional  case  [7],  that  <v>  satisfies  the 
differential  equation : 

[3t  +  (c+e^ic)  3xJ  <v>  =  4  Be2  92<v>  +  i<f>(t,x)<v>  , 


(2.22) 


=  a(x). 
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We  remark  that  the  equation  (2.3)  can  be  derived  by  the  so-called  method 
of  smooth  perturbation  introduced  by  Keller  [1.2].  In  fact,  higher  moments 
of  v  can  be  treated  in  a  similar  manner. 

However,  in  contrast  with  the  traditional  perturbation  methods  for 
moments,  the  current  method  yields  other  important  statistics  than  just  the 
first  few  moments  of  u.  For  instance,  the  statistics  of  the  wave  amplitude 

(2.23)  | u (t , x) | 2  ~  I(t,x)  =  a2[£(t,x)  -  evw(t)] 

2 

can  be  computed  easily.  Suppose  a  (x)  =  p(|x|)  and  p  is  strictly  increasing. 
Then  the  probability 

P{I(t,x)  >  r  )  =  P{|£(t,x)  -  evw(t)  |  <  p”1  (r ) 

1  r+  2 

=  TTJ  /;_exp{ - 2~}  ds 

(27TVE  t)  V£  t 

where  £*  =  £(t,x)  ±  p  J(r)  and  p  1  is  the  inverse  of  p. 

Thus  the  distribution  function  for  I(t,x)  is  given  by 
F(r  ;t,x)  =  1  -  P{I(t,x)  >  r} 

from  which  all  statistics  of  I  at  any  point  (t,x)  is  computable.  Similarly 
the  joint  distributions  at  two  points  Ffr^.r^;  t^,x^,t2,x^) 

=  P{I(tj,x.)  <  r  i ,  i  (t£  ,x.,)  <  r^},  and  multiple  points  may  be  obtained. 

The  computation  becomes  tedious  but  elementary. 

Since  the  phase  0  given  in  (2.19)  is  a  Wiener  functional,  its 
statistics  are  more  difficult  to  compute.  However,  there  exist  an  extensive 
list  of  literatures  on  this  subject,  (e.g.,  see  the  reference  in  [5,8]). 
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3.  Method  of  Progressing  Waves 

p 

Let  Aj  =  Aj(t,x,a)),  j  =  l,2,...,n,  be  some  given  mxn  matrix-valued 
random  functions  of  telk+and  xelR0  ,  depending  on  a  small  parameter  e>0, 
and  let  B(t,x)  be  a  mxn  matrix-valued,  deterministic  function  on  R+  x  IT. 

We  consider  the  random  linear  hyperbolic  system  in  if  for  the  wave  function 
u  in  if  : 

n 

(3.1)  Lu  =  (3  -  l  A.3.  -  B)u  =  0, 

j  =  1  J  J 

8  9 

where  9  =  ,  9.  =  -z —  .  The  random  solution  u(t,x),  in  which  co  is 

t  dt  J  oXj 

surpressed,  is  subject  to  appropriate  initial-boundary  conditions.  For 
clarity  we  shall  present  the  method  of  progressing  wave  solution  to  (3.1) 
in  three  steps  as  indicated  in  §1.  As  a  matter  of  convenience,  a  vector 
will  sometimes  be  regarded  as  a  column  matrix,  such  as  in  (3.1),  and  the 
indicial  notations  will  also  be  used  as  one  sees  fit. 

3.1.  Sample  Progressing  Wave  Solution 

To  study  the  propagation  of  weak  singularities,  such  as  pulses,  the 
progressing  wave  solution  is  an  effective  technique  in  the  deterministic 
case  [9].  For  a  fixed  weft,  we  may  apply  the  method  of  solution  to  the 
system  as  described  in  [10]  and  [11].  To  make  the  paper  self-contained, 
we  shall  present  a  derivation  based  on  Laxfs  approach  [11]. 

For  the  system  (3.1),  we  seek  an  approximate  solution  u  of  the  form: 

(3.2)  u(t,x,w)  =  g[<j)(t,x,u))]v(t,x,aO  , 

where  g  is  a  scalar,  possibly  generalized  function  of  a  real  or  complex 
variable,  <J)  a  scalar  random  function,  while  v  a  m-vector  valued,  random 
function  on  R+  *Rn. 
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Apply  L  to  u  as  defined  in  (3.1)  and  (3.2)  to  get: 


(3.3)  LG  =  g'Mv  +  g Lv, 
where  g'(t)  =  — ,  and 

n 

(3.4)  Mv  =  (<f>  I  -  l  <(>.A.)v  , 

j=l  J  J 

with  (|)^  =  ,  4>_.  =  and  I  a  n*n  identity  matrix.  Since  g  may  be  a 

generalized  function,  for  the  function  Lu  to  be  less  singular,  we  set 

(3.5)  Mv  =  0, 

so  that  (3.3)  becomes 

(3.6)  Lu  =  g Lv  . 

For  (3.5)  to  have  a  nontrivial  solution,  we  must  have 

(3.7)  3t<J>  =  Ae(t,x,3x<[>,a)) , 
and 


(3.8)  v  =  a(t,x,to)r (t ,x,o))  , 

where  X  is  an  eigenvalue,  with  the  corresponding  right  eigenvector  r, 

n 


of  the  random  matrix  £  (f).A.,  and  a  is  a  scalar  random  function.  To  deter- 

j=i J  J 

mine  a,  we  seek  the  next  order  of  approximation  by  writing 

(3.9)  u  =  u  +  q. 

In  view  of  (3.1)  and  (3.6),  the  remainder  q  satisfies 

(3.10)  Lq  =  -g Lv  . 

Assume  the  second-order  approximation  q  to  q  has  a  form  similar  to  u  : 


(3.11)  q  =  h(4>)v1 

where  h  and  are  scalar  and  vector-valued  functions  to  be  determined 
so  that  the  first-order  residual  term  on  the  right-hand  side  of  (3.10)  is 
removed.  This  leads  to  the  choices: 


(3.12)  h'  =  g, 

(3.13)  Mvx  = 


-Lv  . 
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Let  £(t,x,u>)  be  a  left  eigenvector  associated  with  r  so  that  their  inner 
product 

(3.14)  £ • r  =  1  . 

Then,  since  r  is  a  solution  of  the  homogeneous  equation  of  (3.13),  the 
solvability  condition  implies  that 

(3.15)  £ • Lv  =  0 
or,  noting  (3.8) 

n 

(3.16)  [3  -  l  (£-A.r)3.]a  +  (£-L r)a  =  0  , 

z  j=l  33 


which  determines  the  amplitude  function  in  the  first-order  approximation. 

To  find  the  phase  function  <p ,  we  can  integrate  the  Hamilton-Jacobi  equation 
(3.7)  by  the  method  of  characteristics  [10].  The  characteristic  (or  bi¬ 
characteristic)  equations  read: 

'  S  =  -Xq(S-y-q-“>  • 

(3.17)  I,  Hs  ‘  Xyts,y’q'")  ’ 
y(t)  =  x, 


0<s<t , 


^  q(t)  =  . 

In  the  physical  space,  the  characteristic  curve 

(3.18)  ;  y  =  £s(t,x)  ,  o<s<t 

is  called  a  ray  through  the  point  (t,x).  As  in  the  theory  of  geometric 
optics  [12]  ,  the  equation  (3.7)  will  be  called  the  eikonal  equation  for 
the  phase  0,  the  equation  (3.16)  the  transport  equation  for  the  ampli¬ 
tude  a.  It  is  easy  to  check  that,  along  a  ray  F  ,  0  is  conserved  and 

t ,  x 

that  the  transport  equation  (3.16)  can  be  integrated  similarly  to  the  problem 
in  §2,  where  the  characteristic  equation  now  becomes  a  bicharacteristic 
equation  (3.17). 
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3.2.  Diffusion  Approximation 

To  study  the  effect  of  weak  fluctuations,  we  assume  that  the  random 
matrices  in  (3.17)  take  the  special  form: 


(3.19)  Aj  =  A?(t,x,u>)  =  Aj  +  eRj  (t ,x,uj) ,  j=l,2,...,n 


where  A^  is  a  constant  matrix,  and  Rj  a  spatially  homogeneous  random, 
matrix-valued  function  satisfying 


ER.  =  0  , 
i 

where  {p^  are  the  entries  of  R^ ,  i, j ,£,q=l ,2, . . . ,n,  k,p=l , 2 , . . . ,m. 

To  stress  the  dependence  on  e,  the  associated  quantities  will  be  indexed 
by  e,  such  as  £  ,  r  ,  M  ,  etc.  In  view  of  (3.4)  and  (3.14),  we  can  write 


(3.21)  Mere  =  (A£I  -  l  p.A?)rE  =  0 

j=l  J  J 


(3.22)  £eTe  =  1. 


from  which  we  get 


£e-Mere  =  A£  -  £e-(p-A£)r£  =  0  , 


or 


(3.23)  Ae  =  £e- (p-Ae)re  , 


where  p-Afc  =  £  p.Ae  . 


0 

Expand  A  in  the  small  parameter  and  write 


(3.24)  Ae  =  A  +  eA  +  0(e2). 

For  simplicity,  we  put  A  =  A*\  £  =  iP ,  and  so  on.  For  e  =  0  and  by 
(3.19),  (3.22)  and  (3.23)  yield 

(3.25)  £.r  =  1 
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and 

n 

(3.26)  \  =  it '  (P’A)r  =  l  •  (  j>  p.A.Jr. 

j=l  J  J 

By  differentiating  (3.22)  and  (2.23)  with  respect  to  e  at  e=0,  we  get 

n 

(3.27)  A  =  A-(p-R)r  =  *•(  £  p.RJr. 

j=l  J  3 

A 

Note  that  since,  by  assumption  A  is  constant,  the  eigenvectors  l  and 
r  may  depend  only  on  p,  not  on  x.  Similarly,  by  differentiating  (3.22)  and 
(3.23)  with  respect  to  p  and  then  to  x, and  invoking  (3.21)  and  (3 . 24) -  (3 . 27) , 
we  can  express  the  characteristic  equations  (3.17),  keeping  only  terms  up 
to  0(e),  in  the  explicit  form: 


where 


(3.28)  • 


=  Y(q)  +  eY(s,y,q,w) 
=  eQ(s  >y  >q>w)  > 

y(t)  =  x, 

q(t)  =  p. 


(3.29)  ■ 


-  Y.(p)  =  -A'A.r  =  -  I=iAiJk^(p)rk(p), 

n 

Yi(t,x,p,w)  =  -£-Rir  =  -  l  Pif  jk(t,x,u)Jlj  (p)rfc(p) , 

j  7  ^=1 

n 

Qi(t,x,p,w)  =  V  (p*3iR)r  =  l  a^p,  hk(t,x,w)p,«.,(p)L  (p)  , 

j ,k,h=l  J’  J 


i=l , 2, . . . ,n. 

From  the  above  expressions,  we  see  that,  in  general,  the  characteristic 
equations  are  coupled.  In  case  that  they  do  not,  the  asymptotic  analysis 
becomes  much  simpler.  This  special  case  will  be  treated  first. 

(a)  Decoupled  Bi-characteristics 

£  A 

It  is  possible  that  hence,  Y  and  Y  in  (3.28)  may  be  independent 
of  p,  (see  an  example  in  §4).  Then  we  can  decouple  the  first  equation 
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from  the  second  in  (3.28)  to  get  the  physical  characteristic  curve: 


(3.31)  :  £}£  =  Y  +  eY(s,y,ui) ,  Oss<t, 

V(t)=x, 

A 

where  Y  is  a  constant  vector.  This  is  a  special  form  of  the  stochastic 

equation  (1)  treated  in  the  appendix.  Thus  if  the  random  field  Y  (or  R) 

2 

satisfies  the  conditions  of  Theorem  A.l,  then,  as  elO,  t>s+°°  with  T=e  t 
2 

and  o=e  s  held  fixed,  the  solution  of  (3.31)  approaches  the  asymptotic  form 

(3.22)  y(s)  =  5g(t,x)  ~  Cs(t,x)  =  SJt.x)  +  C0CO, 

where 


(3.33)  5s(t,x)  =  x  -  A(t-s), 

and  the  process  {5  (t),  0<o<t}  is  a  diffusion  process  satisfying  the  Ito 
equation  (6) .  Here  the  drift  vector  K  and  the  diffusion  matrix  g,  referring 
to  (A. 6),  are: 


(3.34) 

(3.35) 


Sj,x-AspYj  (T-S2,x-As2)dSjds 

A  A 

,x-As1)Y^ (T-s2,x-As2)ds1ds2 


2* 


i,  j  =  l,2, . . . ,n, 

where  Y^  is  defined  in  (3.29),  and,  by  the  assumption  of  homogeneity, 
K  and  g  are  constant. 


(b)  Coupled  Bi-characteristics 

In  the  general  situation,  Theorem  A.l  does  not  apply  directly  to  the 
bi-characteristic  equations  (3.28).  However,  since  (q-p)=0(e),  if  we 

A 

linearize  Y  about  p, 

(3.36)  Y  (p+q ' )  =  Y(p)  +  Y^'  +  0(e2),  with  q'=q-p, 
then  the  system  (3.28)  becomes 
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(3.37) 


'  g  =  Y(p)  +  Yjq*  +  eY(s ,y ,p+q *  ,uj) 
=  eQ(s,y,p+q’  ,<*>) 


y(t)  =  x, 

■  q’(t)  =  0, 

where 

Yj  =  9qY(q) | q=p 

is  a  nxn  matrix. 

Let 


y*  =  (y.q1). 


(3.38) 


Y* (y*)  =  A*y*  +  C*, 

A*y*  =  (Yjq',0),  C*  =  (Y (p) , 0) , 
and 

■  Y*(s,y*,aj)  =  (Y(s,y,p+q' ,w),  Q(s ,y ,p+q '  ,w) ) . 


Then  the  system  (3.36)  can  be  written  as 


(3.39) 


/  57“  =  Y*(y*}  +  eY*(s>y*>w)> 

y*(t)  =  x*  =  (x,o) 


which  is  of  the  same  form  as  the  standard  equation  (1)  in  the  appendix 
with  y  replaced  by  y*.  Hence  we  can  apply  Theorem  A.l  to  (3.39), 
provided  that  the  corresponding  conditions  (A.l) -(A. 6)  are  met.  To  this 

end  let  K*  be  a  fundamental  matrix  for  the  system  (3.40)  at  e=0  with 

* 

0 


Ki  =  I*.  It  is  simple  to  check  that 


r1  tYii  -i 

Kt  =  [o  I  J  with  <KP  =  K*t 


(3.40) 


Similar  to  Z  as  in  (6) ,  define 


-1, 


(3.41)  Z*(s,z*,u))  =  (K*_s)  Y*[t-s,^*_s(t,x)  +  K*_sz*,aj], 
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where 

(3.42)  q(t,x)  =  (x+Y (p)  (t-s)  ,0)  . 

Then,  according  to  Theorem  A.l,  the  solution  of  (3.40)  has  the  following 
asymptotic  form 


(3.43) 


£*(t,x,e)  ~  £*(t,x)  +  K*_s?*(t), 


where  (£*(t),  0<o<t}  is  a  diffusion  process  in  ]R  4-11  with  the  drift  vector 

and  the  diffusion  matrix  given,  respectively,  by 

t„+T  s,  2nr  ^ 


fAi 

0  J=1 


(3.44) 


T**oo 


'0 


(Z*) . (S, ,C) 


as .v  r ±v  is 

•  j 


(Z*).(s2,C)dSlds2 


VT  fVT 

:o  to 


L  ftf  £(zpi(s1,V(zpj(s2,t)ds14s2. 


The  equation  (3.43)  defines  an  asymptotic  bi-characteristic  curve  in  the 

phase  space.  Its  projection  onto  the  physical  space  yields  the  characteristic 
r,e 

curve  r  : 
t,x 

(3.45)  y  =  CgCt.x)  ~  Cs(t,x)  +  C0(t)  -  (t-sJY^  (t)  =  ^(t,x), 
where 

Cs(t,x)  =  Cs(t,x)  =  x+Y(p) (t-s) , 

(3.46)  { 

1  C*  =  (ca,<£). 

From  (3.45)  we  see  that,  in  contrast  with  the  one-dimensional  case,  the 
random  fluctuation  about  the  unperturbed  characteristic  is  a  diffusion  in 
the  phase,  instead  of  physical  space.  However,  as  to  be  shown  later  in  4§ 
in  a  special  case,  the  physical  characteristics,  or  £  ,  is  asymptotically 
decoupled  from  £ 1  under  appropriate  conditions. 
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3.3.  Asymptotic  Statistics 

In  the  first-order  approximation,  noting  (3.2)  and  (3.8),  the  pro¬ 
gressing  wave  solution  of  (3.1)  takes  the  form 

(3.47)  u(t,x,uj)  ~  aE(t,x,(jo)g[<J>£(t,x,(jd)]r(x) 


where,  within  0(e),  the  unperturbed  initial  eigenvector  r(x)  is  used  in  place 

£ 

of  r  (t,x,Go).  That  is,  r(x)  satisfies  the  equation  (3.21)  with  e=0  and 
Pj=9j0(x).  Suppose  that  the  initial  state  of  u  may  be  represented  as: 

(3.48)  u (0 , x , co )  ~  a(x)g [6 (x)  ]  r(x) , 


with 


(3.49) 


£ 

f  a  (0,t,x,(jj)  =  a(x) , 
-  <J>e(0,x,(jd)  =  0(x) , 

■  rE(0,x,(jd)  =  r (x) . 


By  integrating  the  eikonal  equation  (3.7)  and  the  transport  equation  (3.16) 

£ 

along  the  physical  characteristic  curve  x>  (3.47)  and  (3.48)  yield, 
within  an  error  of  0(e), 


(3.50) 


u (t , x)  ~  a[^E(t,x)]g{0[^Q(t,x)] }r(x) 
X  exp{/J  b[s,£E(t,x)]ds} 


where 

n 

(3.51)  b(t,x)  =  Jt(x) -B(t,x)r(x)  +  £  £  (x)  •  A. 3  .r (x) . 

j  =  l  J  J 

Therefore,  if  the  functions  a,b,g,0  are  such  that  the  conditions  for 
Theorem  A. 2  are  satisfied,  then,  in  the  diffusion  limit,  we  have 


u ( t , x)  ~  a[(o(t,x)]g{0[Co(t,x)])r(x) 
X  exp{/p  b[s,Cs(t,x)]ds}, 


(3.52) 
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where,  depending  on  the  bi -characteristics  being  uncoupled  (case  a)  or 
otherwise  (case  b) ,  the  asymptotic  expression  is  given  by  (3.32)  or 
(3.45),  respectively. 

Suppose  that  (r^  r2>  ...,  rR}  and  (2^,  %v  ....  the  right 

and  left  eigenvectors  correspond  to  the  distinct,  nonzero  eigenvalues 
{\v  A?,  A^},  respectively,  are  such  that 

(3.53)  ^i*rj  =  ^ij *  for  =  l>2,...,k  <  n. 

If  the  initial  state,  instead  of  (3.48),  can  be  represented  as 

m 

(3.54)  u(0,x,u>)  ~  l  a.  (x)g[0  .  (x)]r.  (x) , 

j=l  J  J  J 

then,  by  superposition,  the  corresponding  progressing  wave  solution  (3.51) 
becomes 


(3.55) 


u(t,x,oj)  ~  l  aj[Cjj0(t,x,(i))]g{ej[^j^0(t,x,a))]}rj  (x) 


X  exp{/^  b(s,C j  >g(t  ,x,0)) ] ds} 


where  y=£ .  is  the  asymptotic  solution  of  the  characteristic  equation 
J  ^ 

when 

For  each  component  (mode)  of  u  in  (3.55),  the  asymptotic  statistics  of 
the  progressing  wave  may  be  computed  in  a  similar  manner  as  illustrated  for 
the  model  problem  in  §2.  However,  in  order  to  compute  the  statistics  in¬ 
volving  more  than  one  component  in  (3.55),  one  has  to  obtain  the  joint 
probability  distributions  of  £.  .  This  can  be  done  by  augmenting 

1  ,  S  K  ,  S 

the  single  system  (3.17)  to  k-simultaneous  systems  of  bi-characteristic 
equations,  corresponding  to  the  k  distinct  eigenvalues.  Then  apply  the 
usual  asymptotic  analysis  to  the  simultaneous  systems  to  derive  the  joint 
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diffusion  processes  {£,  , ...,£,  }.  This  involves  no  new  concept,  but 

does  increase  the  computational  complexity. 

We  remark  that,  as  in  the  deterministic  case,  it  is  possible  to 
carry  out  the  higher  order  approximation.  But  the  statistical  problem 
would  become  too  complicated  to  be  of  practical  value.  Also  we  note  that 
the  asymptotic  symbol,  such  as  in  (3.55),  has  a  double  meaning.  That  is, 
the  asymptotic  sequence  is  ordered  by  the  strength  of  singularity  of  the 
phase  function  g  and  the  small  parameter  e .  In  what  follows  we  shall  apply 
the  method  to  some  problems  in  acoustic  wave  propagation. 


4 .  Acoustic  Waves  in  Turbulent  Media 

Consider  the  acoustic  wave  propagation  through  a  weakly  turbulent 
fluid  governed  by  the  following  system 

/  V  *  =  °> 

(4'1)  t  „2 

3tP  +  pc  3^.v  =  0,  t>0,  xelR  ,  n<3, 

where  v  is  the  acoustic  velocity,  p  the  acoustic  pressure,  p  the  density  and 
c  is  the  local  speed  of  sound.  In  this  section,  the  symbol  p  will  be  reserved 
exclusively  for  the  pressure.  The  symbols  3  and  3  •  denote  the  gradient 
and  divergent  operators,  respectively.  Due  to  the  weak  turbulence,  the 
density  p  and  the  sound  speed  c  fluctuate  randomly.  So  we  assume  that  the 
random  functions 


(4.2) 


^  0 

r  P  =  P  (t,x,w)  =  p  +  ep(t,x,w) 

'  c  =  C£(t,X,0))  =  C  +  ec(t,X,(jd) 


depend  on  a  small  parameter 


-  (p V1 

pp  -  Pe(c£)2 


For  convenience,  we  define 


(4.3) 
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so  that 


......  £  £  c  £<2 

(4.4)  =  (c  )  . 


To  put  the  system  (4.1)  in  the  standard  form  (3.1),  let  us  define  u  and 


A.  as  follows: 
J 


(4.5)  u  =  (v1,v2, . . . ,vn,p) , 


(4.6) 


(-D 


0 

0 


0 

0 


£* 

“2Slj 


.  .  .  o 


£  * 

M262j 


Vij 

H62j 


£r 

y.o  . 
1  nj 


\i%6  .  0 

2  nj 


>j  1,2,. . . , n , 


where  6^  =  0  or  1  according  to  i=j  or  i*j 

(4.7)  =  0i  +  ey.  +  0(e2)  ,  i=l,2, 


Let  us  rewrite  (4.3)  as 


where,  in  view  of  (4.2), 


e  A  n .  ^  ^-1  ^  ^  2 

(4.8)  Pj  =  p  ,  y2  =  P (c) 

and 


(4.9) 


\i1  =  -p(t,x,o))/p  , 


y2  =  c[cp(t,x,(jj)  +  2pc(t,x,a))] 


Statistically,  we  assume  that  the  random  functions  p  and  c  are  centered, 
homogeneous  and  satisfy  the  strong  mixing  and  other  smooth  properties  so 
that  the  diffusion  approximation  holds.  For  any  homogeneous  random  functions 
£,n>  let  us  introduce  the  following  notations: 
y  (t,s;x-y)  =  E£(t,x)n(s ,y)  , 

Y|n(t,s;x-y)  =  y^(s,t;y-x)  =  Yn£(t,s;x-y) ,  by  definition. 


(4.10) 
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Then,  by  assumption,  we  have 

(4.11)  Ep  =  Ec  =  0, 
and 

Ypp(t,s;x-y)  =  E{p(t,x)p(s,y)}. 

(4.12)  Ycc(t,s;x-y)  =  E{c(t,x)c(s,y) } , 

Ypc(t,s;x-y)  =  E{p(t,x)c(s,y)}. 

£ 

Now,  noting  (4.7),  the  matrices  Aj  in  (4.6)  takes  the  following  form: 

(4.13)  A*r  =  Aj  +  eR.  +  0(e2)  , 
where 


(4.14) 


Aj  =  A. (y) , 

Rj  =  Aj (y) ,  j=l,2, . . . ,n. 


After  taking  the  above  preliminary  steps,  we  are  ready  to  apply  the 
method  of  §3  to  construct  a  progressing  wave  solution  in  one  and  three 
dimensions . 

(a)  One-Dimensional  Problem 

In  one  space-dimension  (3.5)  and  (3.6)  reduce  to 
(4.15)  u  =  (v,p) 


(4.16)  Ae  =  A(ye)  =  (-1) 


and  (4.13)  becomes 


0  V] 

4  o 


(4.17)  Ae  =  A  +  eR  +  0(e2) 


where,  noting  (4.14) 
A  =  A(y)  = 


(4.18) 


R  =  A(y)  = 


0  V, 
P2  0 

A  , 

0  Vj 

y2  oJ 
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For  the  eigenvalue  problem 


(4.19)  M£(q)r£  =  (y£I  -  qA£)r£  =  0, 


the  eigenvalues  and  the  associated  normalized  right  and  left  eigenvectors  are 


(4.20)  A  £  =  A£j2  =  ±c£q  , 


(4.2!)  r£  =  r£>2  =  ± 


"(p  c  )q 


o  xlr  £  £/, 

»  =  *1,2  =  *2^  C  q  1] 


£  £ 

From  (4.20),  A^  =  ±c  ,  so  that  we  have  the  decoupled  characteristic  equation: 


(4.22) 


&  =  ±c£ 


ds 


±c  (s,y,w)  ,  o<s<t, 


y(t)  =  x  , 


where,  as  is  well  known,  the  different  signs  correspond  to  two  waves 
travelling  in  the  opposite  directions. 

For  simplicity,  let  us  first  consider  only  one  mode,  say  Ae  =  A£  =  ceq . 

A 

To  apply  the  results  in  (case  a)  of  §3,  we  note  that,  by  (4.2),  Y  and  Y  in 
(3.31)  are  given  simply  by 

Y  =  c  , 

(4.23) 

Y  =  c(t,x,oj)  . 

Thus  the  asymptotic  ray  (3.32)  becomes 

(4.24)  y  =  £s(t,x)  =  x  -  c(t-s)  +  ^(t)  , 


where  satisfies  the  Ito  equation 


(4.25) 


dc  =  <da  -  vdw(x-a) ,  with  3=v  , 


cT  =  °  • 


The  parameters  k  and  3  are  defined  as  in  (3.34)  and  (3.35) 

t^+T  s 

1v»  fT-c  T-c  -c  'llrU  rls7 


k  =  lim  ^  /  °  /  1 y'  [T-S  ,T-s  ;c(s  -s,)]ds  ds. 

T-hx  1  t0  t0  cc  1  1  Y  1  1  ‘ 

t  +T  t  +T 

3  =  lim  y  lt°  ftU  Ycc[T-s1,T-s2;S(s1-s2)]ds1ds2 


T-»oo 


(4.26) 
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in  which  y^c(t,s;x)  =  9xYcc(t,s;x) . 


Let 

the  initial  condition  be  given  by 

(4.27) 

u(0,x,w)  ~  a(x)g[0(x)]r(x) . 

A 

Here  r(x)  is  the  initial  right  eigenvector  of  (0^A) .  As  seen  from  (4.21), 


for  e=0 , 

the  resulting  eigenvectors 

(4.28) 

r  - l*i 

-(pc)  , 

„  1»-  aa  _ 

r  =  ,  A  =  j  [-pc  1] 

1  J 

are  constant.  Consequently,  the  function  b  defined  by  (3.51)  vanishes,  or 
(4.29)  b(t,x)  E  0. 

Hence  the  asymptotic  solution  (3.52)  reduces  to 


(4.30) 

r~  aa  —  1  i 

-(Pc) 

u(t,x,w)  ~  a[Co(t,x)]g{0[5o(t,x)]} 

from  which  the  second  component  gives  the  acoustic  pressure 

(4.31)  p(t,x,w)  ~  a[£ (t,x) ]g{0 [£ (t,x) ] } 

where  £(t,x)  =  £^(t,x)  is  the  diffusion  process  given  in  terms  of  the 
Wiener  process  {w(t),t>0}  : 

(4.32)  £(t,x)  =  x  -  (c+e2K)t  +  evw(t) . 

Now  it  is  clear  that,  if,  instead  of  (4.27),  the  initial  state 
consists  of  two  modes: 


(4.33) 

2 

u(0,x,w)  ~  l  a. (x)g. [0. (x)]r. 

j=l  J  J  J  J 

the  pressure  field  is  simply 

2 


(4.34) 

p(t,x,w)  ~  l  a. [£, (t,x)]g.{0.  [£. (t,x)] }, 
j=i  J  J  J  J  J 

where 

(4.35) 

Cj(t,x)  =  x  +  (-1)-^  [(c+e2K)t  +  evw(t)]. 
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Therefore  the  statistics  of  p  are  certain  Wiener  functionals  which  can 
be  evaluated  in  exactly  the  same  manner  as  that  of  the  model  problem  in 
§2.  Also  it  is  interesting  to  observe  that  the  asymptotic  pressure  wave  p 
given  by  (4.30)  depends  only  on  the  fluctuation  of  the  local  sound  speed 
c  ,  but  not  on  the  density  fluctuation  p,  as  assumed  in  (4.2).  This  also 
turns  out  to  be  true  in  higher  dimensions  as  to  be  shown  below. 


(b)  Three-Dimensional  Problem 
In  three  dimensions,  we  set 


(4.36) 


(4.37) 


u  =  (vx >v2,v3,p) 


Aj  =  A.(y)  =  (-1) 


0 

0 

0 


0 

0 

0 


0 

0 

0 


£  r  £  O  £  ^ 

LM2S1 j  *2S2J  *2«3J 


uiau 

£x 

^l62j 

v\. 

0 


Analogous  to  (4.20)  and  (4.21),  the  non-zero  eigenvalues  and  the  eigen- 
3 

vectors  for  (  £  q.A.)  are  found  to  be 

j-1  3  3 


(4.38) 


\e  =  \C 


1,2 


±ct|q| 


(4.39)  r 


where 


/v  /  e  e 
-q/p  c 


1,  e  e*  , , 

J[-P  C  q  1]  , 


(4.40)  q  =  q/lql=  -|— y  (q1,q2»q3)  . 
From  (4.38),  we  have 


(4.41)  a/  =  aq/2  =  , 

C4.42)  -  »/>2  .  41,13/  . 
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Thus  the  bi-characteristic  equations  (3.28)  takes  the  form: 


(4.43) 


dy  rAA  .  \  A  1 

=  +[cr  +  ec(s,y,w)r]  , 


=  ±e  3y.c(s,y,aj)  |r|  ,  0<s<t, 

y(t)  =  x, 
r(t)  =  q. 


which  is  of  the  case  (b)  in  §3.  To  apply  the  results  obtained  there,  let 
(4.44)  r  =  q  +  q' . 

Then  by  linearizing  the  term  cr  in  (4.43),  we  get  the  system  (3.37),  where 


(4.45)  J 


Y(q)  =  +cq  , 

V* =  +  (q,,;)"}  ’+ 1 

Y(s,y,q+q' ,w)  =  +c(s,y,w)|3lH_|.  , 
Q(s,y,q+q' ,w)  =  ±9yc(s,y,u>)  |q+q' |  . 


Introduce 


(4.46)  y*  =  (y,q ' )  =  (yj,...,yg)  , 


and  define  Y*,  A*  and  Y*  as  in  (3.38),  and  Z*  as  in  (3.41).  That  is 


(4.47)  Y 


Y.  =  4C(s,y,w) |q+q' |"1(q+q') . ,  j=l,2,3. 


J  *  j  *  +cls ,y ,u)j | q 

=  ±9^._3c(s,y. 


w)|q+q'|  ,  j=4,5,6. 


(4.48)  (Z 


*).  =  P 

tb  lv, 


j(t-s,-)  -  (t-s)Yj  +j(t-s, • )  , 


Yj (t-s ,  •  )  , 


j=l,2,3  , 
j=4,5,6. 


where 


(4.49)  Yn.(or  Y!)(t-s,-)  =  Y.  (or  Y! )  (t-s  ,SQ  (s  ,x)  +  y  +  (t-s)q' ,q' ,oi)  , 


J  '  J 


J  v  3 


Cs(t,x)  =  x  ±  cq(t-s) . 


(4.50) 
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To  compute  the  drift  vector  and  the  diffusion  matrix  defined  by  (3.44), 
we  make  use  of  (4.45) -(4.50) .  In  order  to  express  the  results  concisely, 
let  us  introduce  the  function  bij  as  follows: 


(i)  For  l<i,j<3. 


✓V 


A., 


bij^Sl,S2,q+ql)  =  Ycc^S1,S2;c(q+q,) (S2"S1^ (q+q ' ) i Cq+q ' ) j 


(4.51) 


*  tsrs2)3iYcc 11  ‘‘'♦‘''b 


-  s 


l's23i3jYcc  l> 


[si »s_;y]  at 


where  we  set  S.y  [•]  =  7.  , 

i'ccL  J  3y.  'ccL  1  2 


y  =  c(q+q') (s2-s1) 


(ii)  For  l<i<3,4^j<6, 

(4.52)  (s1,s2;q+q')  =  3j _3ycc N (q+q ' ) t  +  s18i3jYcc l  •  ] I q+q ' I 

(iii)  For  4<i,j<6, 

(4.53)  bij  (s1,s2Jq+q')  =  "3i-38j-3Ycc[‘] iq+q' I2  • 

In  terms  of  b^j  ,  the  equations (3. 44)  yield: 

1  tn+^  tn+^ 

6ii (q+q*)  =  lira  f  !t  jt  b ., (s  ,s  ;q+q')ds,ds 
J  t-h»  0  0  J 

t.+T  s,  3 


(4 


54)  Ki(q+q')  =  lim  i  /  °  //  {  [  3.b  (s,,s  ;q+q') 

T-*»  0  0  j  =  l  3  J  z 


JA  3^rTbij(si’S2;q+q,:i}dslds2’  1  >  j  =1  >  2  > 

J=4  nj-3 


,6. 


Assume  all  of  the  above  limits  exist.  Then  for  large  s^  and  s„ ,  in  view 
of  (4.51)  -  (4.53),  the  components  b^j  in  (4.51)  with  lsi,j<3  are  dominani 
If  the  corresponding  limits  are  finite,  the  rest  must  be  zeros.  That  is. 


(4.55)  =  0,  6i j  =  0  for  i>4  or  j>4. 
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This  has  an  interesting  consequence.  Namely,  in  the  asymptotic  approxi¬ 
mation,  the  process  {y*(s),  0<s<t}  becomes  degenerate  and  the  time-averaging 
in  (4.54)  has  the  effect  of  projection 


(4.56)  y*(s)  (y(s)  ,0) 


from  the  phase-space  into  the  physical  space,  where 


(4.57)  y(s)  ~  5s(t,s)  =  Ss(t,x)  +  ^(t). 


Here  £  is  a  diffusion  in  fl3  defined  by  the  Ito  equation 


(4.58) 


,  dr  .  =  K.da  +  7  V- -dw. 

a, i  i  iJ  J 

■  «t.i  ■  0  ■ 


(t-o) , 


i=l , 2,3 


3 

in  which  K.  =  K.(q)  and  V  V-.V.  •  =  8,,  •  (q)  with  k  defined  by  (4.54) 

1  1  =  IK  KJ  KJ  1  1J 

for  i,j=l,2,3. 


Suppose  the  initial  state  has  only  a  single  mode  corresponding  to  the 
first  eigenvalue  A^  =  ce|q|. 


(4.59)  u(0,x,w)  ~  a(x)g[0(x)]r(x) , 

where,  noting  (4.21) 


(4.60)  r (x)  =  rx(x) 


A  A  _  1  A 

-(pc)  q 


J  > 


is  the  associated  unperturbed  right  eigenvector  with  the  left  counterpart 
given  by 


(4.61) 


*(x)  =  t  [-PCR  i] »  q  = 


 ve (x) 


rvei 


With  the  aid  of  (4.60)  and  (4.61),  we  can  easily  compute  b  defined  by 
(3.51)  to  get 

1  3 

(4.62)  b(x)  =  ■=■  l  I  3.9.  (x),  with  0.  =  3.0/1701, 

j=l  J  J  J  J 
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which,  in  contrast,  is  zero  in  one-dimension.  This  term  reflects  the 
curvature  effect  of  the  initial  wave  front.  In  view  of  (4.57),  (4.58)  and 
(4.62),  the  asymptotic  progressing  wave  satisfying  the  condition  (4.59)  has 
exactly  the  same  form  as  that  of  (3.52) 


u(t,x,to) 

(4.63) 


a[£0(t,x)]g{e[£0(t,x)] } 


V0  (x) 


X  exp  {/q  b[£s(t,x)]ds}  , 


from  which  we  get  the  acoustic  pressure  field: 

(4.64)  P (t , x)  ~  a[£0(t,x)]g{e[£0(t,x)]}  exp  {/*  b[£s(t,x)]ds}. 


It  is  instructive  to  work  out  two  concrete  examples  for  which  the  phase 
factor  g  has  a  jump  discontinuity: 

fl,  t>0 

(4.65)  gCt)  =  H(t)  = 

^0,  t<0  , 

where  H  is  the  Heaveside  Function.  Let  us  consider  the  following  special 
cases  of  the  phase  function  0(x)  corresponding  to  the  plane  and  the  spherical 
waves ,  respectively. 


(case  i)  The  Plane  Wave 

In  this  case,  we  assume 


(4.66)  6  (x)  =  orx,  |a|  =  1. 

The  above  implies  that  q  =  V0  =  a  and  b  =  0. 

The  pressure  wave  (4.64)  reduces  to 

(4.67)  P ( t , x)  ~  a[£0(t,x)]H[a*£0(t,x)], 
where 

£0(t,x)  =  X  -  C£at  +  eNw(t) , 

(4.68)  2 

c  =  c  +  e  K(a),  N  =  [v. -r  with  NN*  =  6(a), 
t'  1 J 

3 

and  w(t)  is  the  standard  Wiener  process  in  E  . 
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The  drift  vector  k  and  the  diffusion  matrix  g  are  determined  by  (4.54). 

The  statistics  of  p  can  be  computed  easily.  For  instance,  let 

(4.69)  I (t , x)  =  | p (t ,x) 1 2 

be  the  intensity  of  p.  Then,  from  (4.67),  the  mean  and  the  covariance  of  I 
can  be  shown  to  be: 


I(t,x)  =  E{ I (t ,x) } 


(4.70) 


“5?  W.x)|a(x-V‘t*Nz)|2  exP{TF-)dz. 


✓v  A 


and,  for  0<t<s, 

(4.71)  Cov.{l(t,x),  I (s ,y) }  =  M4(t,s;x,y)  -  I (t,x) I (s,y) , 
where 


M4(t,s;x,y)  =  E{ I (t ,x) I (s ,y) } 

(4’72)  1  ,  ,  .  .  |2 

=  - 3 - FT  'Bft  x)  'B(s  yJa(x-c  at+Nz)| 

(2ttj[t(o-t)]  ^  1 

9  ^  I  i  I  ^ 

X  Ia(y-C£as+Nz')  |  exp{- 
Here  B(t,x)  denotes  the  half-space: 

(4.73)  B(t,x)  =  {zeIR3:  a  *  (x-c^at+Nz)  >  0}. 

In  fact  the  multi-point  statistical  averages,  such  as 

$^t1»t2,’*’;x^1hx^2h--.)  =  E{F[p(t1  ,x^h  ,p(t2,x^-) , .  .  .] }  , 

may  be  computed  just  as  easily,  since  we  know  all  the  finite-dimensional 
joint  probability  distributions  of  the  Wiener  process. 

(case  ii)  The  Spherical  Wave 
Suppose 


(4.74)  0 (x)  =  r 


0 
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Then 

Y  ^  X  ^ 

V0  =  - =■  ,  V0  =  -r— i-  =  x  , 

ur  11 

(4.75)  and 

~  2 
V-(V0)  «  ^  . 

Thus,  by  (4.62) 

(4.76)  b(x)  =  y£y  . 

The  pressure  field  (4.64)  now  becomes 


p(t,x)  ~  a[50(t,x)]H{r0-lC0(t,x) I } 

(4.77)  A  t  . 

x  exp{c  /Jj  ICs(t,x)  I"  ds} 

where 

(4.78)  £s(t,x)  =  x  -  c£x(t-s)  +  N[w(t)-w(s)]  , 

(4.79)  c  =  c  +  e2K(q),  NN*  =  B(q) 
with 


Therefore,  unlike  the  previous  case,  the  drift  <  and  the  diffusion  matrix 
3  are  now  functions  of  x.  This  is  due  to  the  variation  of  the  normal  vector 
field  to  the  wave-front:  6(x)=const.  Also,  due  to  the  wave-front  curvature 
effect  (b*0),  the  average  intensity  I  can  no  longer  be  calculated  explicitly 
as  done  in  (4.70).  Instead,  noting  (4.77),  we  set 


I(t,x) 

(4.80) 


E{a2U0(t,x)]H[r0-|£0(t,x)  |] 

x  exp [2c  /q  Us(t,x)  r1ds]  }  . 


which,  in  view  of  (4.78),  is  an  average  Wiener  functional  or  a  Wiener 
integral.  To  evaluate  (4.80)  approximately,  one  may  introduce  an  asymptotic 
expansion  in  e,  or  adopt  the  method  of  differential  equation,  [13],  [5]. 
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In  the  former  approach,  the  theory  of  large  deviation  [14]  could  be  useful. 
Alternatively,  we  regard  the  average  in  (4.80)  as  taken  with  respect  to 


the  diffusion  process  £s(t,x)  with  the  drift  c^fx)  and  the  diffusion 
matrix  B(x).  By  applying  the  Eeynman-Kac  formula  to  (4.80)  along  the 


effective  characteristic,  one  can  show  that  I  satisfies  the  following 
equation : 


L,Wx)y] ♦ 


l(0,x)  =  |a(x) |2H(r0- |x|)  . 


to  which  one  can  apply  some  method  of  asymptotic  solution  for  parabolic 


equations  to  yield  an  approximate  evaluation  of  I.  It  is  also  possible  to 
derive  differential  equations  for  higher  moments  at  one  time,  such  as 
defined  in  (4.71)  with  t=s.  However  this  will  not  be  done  here.  A  similar 
procedure  was  discussed  in  our  papers  [5,8] . 

Concluding  Remarks 

In  closing  we  wish  to  make  the  following  remarks: 

(i)  The  method  of  progressing  waves  introduced  in  this  paper  has  a  wide 
range  of  applications.  Even  though  we  have  only  discussed  its  appli¬ 
cation  to  acoustic  problems,  it  is  clear  that  the  method  applies  to 
other  problems  in  electromagnetic  waves  and  other  waves  in  continuous 
media. 

(ii)  One  of  the  shortcomings  of  the  method  is  the  linearization  of  the 


mean  vector-field  Y(q)  in  the  bi-characteristic  equations  (3.28)  in 
order  to  apply  the  existing  diffusion  limit  theorems.  But  this 
excludes  the  strong  interaction  between  the  mean  and  fluctuating 
fields.  The  remedy  to  this  problem  requires  a  new  type  of  limit  theorem 
that  generalizes  the  results  of  Kesten  and  Papanicolaou  [15] . 
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(iii)  As  shown  in  the  three-dimensional  acoustic  problem  (b)  in  §4,  the 
linearization  of  Y  leads  an  asymptotic  decoupling  of  the  bi¬ 
characteristic  equations.  This  simplifies  the  subsequent  statistical 
computation  considerably.  In  fact  this  asymptotic  reduction  holds 
true  in  the  general  case  as  presented  in  §3. 

(iv)  In  view  of  the  special  cases  (i)  and  (ii)  in  §4,  the  statistical 
problem  for  the  plane  wave  is  much  simpler  than  the  spherical  one, 
or,  for  that  matter,  any  non-plane  waves.  This  suggests  that  the 
method  of  plane  wave  [16]  may  be  used  for  simplifying  the  statistical 
computation . 

(v)  In  developing  the  progressing  wave  solution,  we  are  mainly  interested 
in  propagation  of  weak  singularity.  For  continuous  waves,  e.g., 

g(0)  =  exp{k0},  it  may  be  interpreted  as  an  asymptotic  expansion  in  the 
reciprocal  powers  of  the  large  parameter  k  [17] .  For  time-independent 
media,  the  procedure  yields  the  well-known  geometric  optics  method. 

The  asymptotic  analysis  of  the  stochastic  ray  equations  near  a 
caustic  was  done  by  White  [18]. 
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Appendix 


Consider  the  random  characteristic  equation  in  R  n: 

[  7T7  =  0<s<t , 

(1)  |  ds 

y(t)  =  x  e  Rn. 

0  0  0 

where  Y  =  (Yj,...,Yn)  is  a  given  family  of  random  fields  for  ee[0,E.g], 
defined  over  Rx  Rn  x  ft.  Let  F  be  o-field  of  measurable  sets  in  the 
sample  space  ft  on  which  a  probability  measure  P  is  prescribed. 

Suppose  that  Y  satisfies  the  conditions 

(A.l)  Ye(t,x,u>)  =  Y(t,x)  +  eY(t,x,u)), 


where 


Y(t,y)  =  A(t)y  +  c(t)  . 

and  A(t)  is  a  nxn  matrix,  and  c(t)eRn  are  bounded  and  continuous.  Let 

A 

y  =  £s(t,x)  be  a  solution  of 


(2)  1 

,  ft  -  ?(s.y) 

1 

k  y(t)  =  x. 

Set 

(3) 

y  =  Ss(t,x)  +  Ks(t)z, 

where  Ks(t)  is  a  fundamental  matrix  for  the  system  (2)  with  c(t)=0  and 
K  (t)=I.  Then,  in  view  of  (1)  and  (2),  it  is  easy  to  verify  that  z(s) 
satisfies 


(4) 


4^-  =  eK  ^[s,^  +K  z,u>]  , 
ds  s  1  ’^s  s  J 

z(t)  =  0. 
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Let  us  set 

(5) 


Zs  =  z(t-s). 


(6)  Zt(s,z,w)  =  K^_s(0)Y[t-s,£t_s(t,x)+Kt_s(0)z,w] . 


Then  (4)  yields  an  initial-value  problem: 


(7) 


dzs 

Jg-  =  eZt(s,zs,w)  ,  0<s<t, 


zo  ’  °- 


In  this  form,  we  can  apply  the  limit  theorem  of  Khasiminskii [3]  or,  more 
generally,  that  of  Papanicolaou  and  Kohler  [19].  To  this  end  we  assume  that 


(A2)  There  exists  a  family  of  a-field  for  Q  3 
^1  ^2 

F  C  Fs2  C  F’  V[Sl5tl]  C  ^S2jt2^  C  K 

(A3)  is  strongly  mixing  in  the  sense  that  g  p(t)  3 

sup.  sup  | p (A |  B)  -  p(A)  |  =  p(t)  4  0  as  tt00  , 

00 
r 

s+t 

7S 
-oo 


seR  A  eF 

c 

BeFS 


and 


Vi 


/oo 

Qp  (t) dt  <  0°. 


(A4)  The  random  field  Y  in  (A.l)  is  jointly  measurable  in  its  arguments 
and,  for  fixed  (t,x),  Y(t,x,  )  is  F^-measurable . 

(A5)  There  is  an  absolute  constant  C>0  3  |Y(t,x,co)|  <  C(l+|x|), 

and  the  partial  derivatives  up  to  order  four  are  bounded,  or 


9aY(t,x,to) 


Bx 


a 


<  C,  for  |a| =1 , . . .,4  , 


where  a=(a^, . . . ,ou)  denotes  a  multi-index. 
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(A6)  EY(t,x)  =  0, 

and  the  following  limits  exist 

.  to+T  to+T  . 

e.,(z)  =  lim  ±  /  J  EZ'(s  ,z)z}(s  z)dsds 
J  0  0 

t_+T  s.  ~„i ,  . 

l  0  In  3ZT(s1>z) 

IC.  (z)  =  lim  =  /  J  E  l  - 7T- -  zj,(s  ,z)ds  ds  i,j=l,...,n 

T-*»  0  0  j  =  l  j 

which  are  independent  of  t  ,  and  (t,x),  (see  (6)).  Here  Z*  denotes  the 
it^1  component  of  Z 

Under  the  above  assumptions,  we  can  easily  check  that  the  assumptions 
of  the  theorem  of  Papanicolaou  and  Kohler  [19]  are  satisfied.  Hence  we 
have 


Theorem  A. 1 .  Let  the  conditions  (A.l)  to  (A6)  be  fulfilled,  and  let 

2  2  e 

a=e  s,  x=e  t.  Then,  for  0<c<t<T,  the  solution  process  zs=£s(t,x)  converges 

weakly  as  e+O  to  a  diffusion  process  C0(T)  with  Ct(t)=0,  for  which  the 

drift  vector  k(z)  and  the  diffusion  matrix  3(z)  are  given  in  (A6)  . 

n 

Remarks:  (i)  Let  us  set  3. . 

-  v  ij 

A 

process  £^(x)  will  satisfy  the  Ito  equation 

(  dc  =  <(C_)d a  +  N (C) dw(x-a) , 

(8)  \  0  °  ° 


1  v.,v,  .  or  B=NN*.  Then  the  limiting 
k=l  lk  kJ 


1  q =  0 

where  N=(v^)  and  w(t)  is  a  standard  Wiener  process  in  E 


(ii)  As  a  corollary,  in  view  of  (3),  (5),  we  have 


y(s)  ~  Cs(t,x)  +  ks;o(t) • 

(iii)  If  A(t)S0,  c(t)  is  independent  of  t  in  (A.l), 

A 

then  £s=x+c (t-s)  ,  Ks=I  (identity). 


(iv)  If  Y(t,y)  in  (A.l)  is  nonlinear,  a  diffusion  limit  may  not  exist. 

For  a  special  case,  one  is  referred  to  the  paper  [19]  by  Kesten  and 
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Papanicolaou  concerning  a  stochastic  acceleration  problem. 

To  apply  the  above  theorem  to  a  stochastic  wave  problem,  we  must 
ensure  the  weak  convergence  of  a  certain  function  or  functional  of  the 
solution  process  y=£^(t,x)  of  the  system  (1),  such  as  the  amplitude  and 
the  phase: 

ab(t,x)  =  a[^(t,x)]  , 

(9) 

eE(t,x)  =  /0b[s,£g(t,x)]ds. 

Suppose  the  following  conditions  hold: 

(A. 7)  There  exist  positive  numbers  ot,(3,M  3 

EICqC^.x)  -  SQ(t2,x)|a  <  Ml t1“t2 1 1+e ,  V  t1,t2>0, 

(A. 8)  lim  lim  sup  P{|££  (t,x)-£C  (t,x)|>6}  =  0 
hlO  e4-0  |  s ^ - s ^  I  ^h  S1  S2 

V(t,x),  6>0. 

(A. 9)  a  :  ]R n  •+•  F  and  b  :  R+*Rn  -*■  Rare  bounded  continuous  and 
3  g  :  1R  n  •+  R+  3 

| b(t,x) | 

llm  SUP  pry-)  =  °- 

N-*»  t>0  81  ' 

|  x  |  >N 

Then  we  have  the  following  asymptotic  result: 

Theorem  A. 2 .  Let  the  conditions  (a. 7) -(A. 9)  be  satisfied.  Then,  in  the 
diffusion  limit  as  in  Theorem  A.l,  we  have,  referring  to  (9), 

(a)  aff^Ct, x)]  ~  a[£g(t,x)]  , 

(b)  0b(t,x)  ~  9(t,x)  =  /ob[s,^s(t,x)]ds, 
where  £®(t,x)  ~  £s(t,x)  =  £s(t,x)+MsCo(t) . 

Remarks :  (i)  Here  the  asymptotic  equality  M~",  instead  of  the  weak 

convergence  means  the  weak  convergence  with  respect  to  the  centered 
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£  ^ 

process  E,^  (t  ,x)  -^s  (t  ,x) ,  since  the  latter  part  may  not  have  a  limit  as 
s<tf°°. 

(ii)  The  continuity  assumption  on  the  function  a  may  be  relaxed  to  the 

so-called  a.s.  B-continuity  (for  definition,  see  p.282,  [20]).  This 
will  allow  a  to  have  simple  jump  discontinuities. 

The  proof  of  the  theorem  follows  from  two  theorems  in  Gikhman  and 
Skorokhod  [21],  (Theorem  2  on  p.  450,  and  Theorem  2  on  p.  486  with  the 
obvious  error  h-*»  in  (3)  corrected  to  h+0) .  The  generalization  from  their 
version  in  one-dimension  to  one  in  several  dimensions  is  immediate. 
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